Microscopic universality in the spectrum of the lattice Dirac operator 
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Large ensembles ol complete spectra of the Euclidean Dirac operator for staggered fermions are 
calculated for SU(2) lattice gauge theory. The accumulation of eigenvalues near zero is analyzed 
as a signal of chiral symmetry breaking and compared with parameter-free predictions from chiral 
random matrix theory. Excellent agreement for the distribution of the smallest eigenvalue and the 
microscopic spectral density is found. This provides direct evidence for the conjecture that these 
quantities are universal functions. 
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Hadronic properties such as the lightness of the pion 
masses and the absence of parity doublets, strongly in- 
dicate that chiral symmetry is broken spontaneously. In 
QCD, a great deal of insight in such nonperturbative phe- 
nomena has been obtained from extensive lattice QCD 
simulations |l|,||. This did not go without a significant 
amount of effort. One of the difficulties is that the order 
parameter of the chiral phase transition, {iptp}, can be ob- 
tained only after a complicated limiting procedure: the 
thermodynamic limit, the chiral limit, and the contin- 
uum limit. In addition, in the chiral limit it is extremely 
costly to take into account the effect of the fermion deter- 
minant. Since the fermion determinant can be expressed 
as a product over the Dirac eigenvalues, this alone war- 
rants a detailed study of the QCD Dirac spectrum. More- 
over, (tpip) is directly related to the QCD Dirac spectrum 
through the Banks-Casher relation 



(ijnp) = hm lim ^(0). 



(1) 



Here, m is the quark mass, V is the volume of space-time, 
and p(A) = (5Z„£(A ~ ^«)) 1S the eigenvalue density of 
the Euclidean Dirac operator, iD = ij^d^ + j^A^, aver- 
aged over gauge field configurations. We observe that the 
average position of the smallest eigenvalues is determined 
by the chiral condensate. In this letter we focus on fluc- 
tuations of the smallest eigenvalues about their average 
position. It should be clear that such fluctuations affect 
the fermion determinant and are important for the un- 
derstanding of finite size effects [Q . The hope is that they 
are given by universal functions which can be obtained 
analytically. This analytical information could then be 
used to facilitate extrapolations to the thermodynamic 
and chiral limits. 

A similar situation arises in mesoscopic physics Q . In 
these studies, it was shown that for a sufficient amount 
of disorder, spectral correlations are universal and can 
be obtained from a Random Matrix Theory (RMT) with 



only the basic symmetries included. On the other hand, 
the average spectral density is non-universal and requires 
specific knowledge of the dynamics of the system. In this 
letter we investigate the question whether a similar sepa- 
ration of scales takes place in QCD. Does the disorder of 
lattice QCD gauge field configurations result in universal 
fluctuations of the small Dirac eigenvalues? 

According to the Banks-Casher relation, the low-lying 
Dirac eigenvalues are spaced as 1/V for (•tptp) ^ 0. Recent 
work by Leutwyler and Smilga || shows that this part of 
the spectrum is related to the pattern of chiral symmetry 
breaking by means of a class of sum rules for the inverse 
Dirac eigenvalues. It is natural to magnify the spectrum 
near A = by a factor of V. This leads to the intro- 
duction of the microscopic spectral density defined by 



p s (z) = lim 

V —>oo 



(2) 



where £ is the absolute value of {ipip). Based on the 
analysis of the Leutwyler-Smilga sum rules, it was con- 
jectured Q that this distribution is universal and only 
determined by the global symmetries of the QCD parti- 
tion function, the number of flavors, and the topological 
charge. If that is the case it can be obtained from a much 
simpler theory with only the global symmetries as input. 
Such a theory is chiral RMT which will be discussed be- 
low. Whether or not QCD is in this universality class is 
a dynamical question that can only be answered by lat- 
tice QCD simulations. The investigation of this question 
is the main purpose of this letter. At this moment it 
can only be addressed on relatively small lattices where 
our results are consistent with zero topological charge. 
The pertinent question of what happens in the contin- 
uum limit has to be postponed to future work. In this 
limit we expect zero modes to become important. p s {z) 
is then different in different sectors of topological charge. 
However, on present day lattices with staggered fermions 
there seems to be no evidence of a "zero- mode zone" M , 
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and the situation is controversial at best. 

There are already several pieces of evidence supporting 
the conjecture that p s is universal: (1) The moments of 
p s generate the Leutwyler-Smilga sum rules 0. (2) p s is 
insensitive to the probability distribution of the random 
matrix 10 (3) Lattice data for the valence quark 
mass dependence of the chiral condensate could be un- 
derstood using the analytical expression for p s fll2| , |l3| ] . 
(4) The functional form of p 8 does not change at finite 
temperature |l4|. (5) The analytical result for p s is found 
in the Hofstadter model for universal conductance fluc- 
tuations | p5| . (6) For an instanton liquid p s shows good 
agreement with the random-matrix result [fl6|| . However, 
a direct demonstration for lattice QCD was missing. 

An analysis of Dirac spectra on the lattice was per- 
formed in Ref. jL7| where it was shown that the spectral 
fluctuations in the bulk of the spectrum on the scale of 
the mean level spacing are universal and described by 
RMT. This showed that the eigenvalues of the Dirac op- 
erator are strongly correlated. Only few configurations 
were available in this study, but spectral ergodicity al- 
lowed to replace the ensemble average by a spectral av- 
erage. However, spectral averaging is not possible for p s 
since only the first few eigenvalues contribute. Therefore, 
a large number of configurations is essential. 

We briefly summarize the main ingredients of chiral 
RMT. In a random-matrix model, the matrix elements of 
the operator under consideration are replaced by the ele- 
ments of a random matrix with suitable symmetry prop- 
erties. Here, the operator is the Euclidean Dirac operator 
iD which is hermitian. Because 75 anti-commutes with 
iD the eigenvalues occur in pairs ±A. In a chiral basis, 
the random-matrix model has the structure Q 

W' 



iD 



im 



where W is a matrix whose entries are independently 
distributed random numbers. In full QCD with Nf fla- 
vors, the weight function used in averaging contains the 
gluonic action in the form exp(— Sgi) and Nf fermion de- 
terminants. In RMT, the gluonic part of the weight func- 
tion is replaced by a Gaussian distribution of the random 
matrix W. The symmetries of W are determined by the 
anti-unitary symmetries of the Dirac operator. Depend- 
ing on the number of colors and the representation of the 
fermions the matrix W can be real, complex, or quater- 
nion real |]l8| . The corresponding random-matrix en- 
sembles are called chiral Gaussian orthogonal (chGOE), 
unitary (chGUE), and symplectic (chGSE) ensemble, re- 
spectively. The microscopic spectral density has been 
computed analytically for all three ensembles ]9|, [l9| , p0t . 

We have performed numerical simulations of lattice 
QCD with staggered fermions and gauge group SU(2) 
for couplings /3 — 4/g 2 — 2.0, 2.2, and 2.4 on lattices of 
size V = L 4 with L = 8, 10, and 16. This range of lat- 
tice parameters covers the crossover region from strong 



to weak coupling of SU(2) pl| . The boundary conditions 
are periodic for the gauge fields and periodic in space and 
anti-periodic in Euclidean time for the fermions. In this 
work, we only study the quenched approximation using a 
hybrid Monte Carlo algorithm [^2) . This made it possible 
to generate a large number of independent configurations 
(indicated in the figures). The analysis of unqucnchcd 
data with 4 dynamical flavors is in progress. 

In SU(2) with staggered fermions, every eigenvalue of 
iD is twofold degenerate due to a global charge conjuga- 
tion symmetry. In addition, the squared Dirac operator 
— D 2 couples only even to even and odd to odd lattice 
sites, respectively. Thus, — D 2 has V/2 distinct eigenval- 
ues. We use the Cullum-Willoughby version of the Lanc- 
zos algorithm |23| to compute the complete eigenvalue 
spectrum of the sparse hermitian matrix — D 2 in order to 
avoid numerical uncertainties for the low-lying eigenval- 
ues. There exists an analytical sum rule, tr(— D 2 ) = 4V, 
for the distinct eigenvalues of — D 2 (^4|. We have checked 
that this sum rule is satisfied by our data, the largest rela- 
tive deviation was ~ 10~ 8 . We have also made a detailed 
study to determine the optimal acceptance rates and tra- 
jectory lengths . The integrated autocorrelation times 
are in the range of 1 to 4. The chiral condensate was ob- 
tained by fitting the spectral density and extracting p(0) 
and is given in Table || below. 

The overall spectral density of the Dirac operator can- 
not be obtained in a random-matrix model since it is 
not a universal function. The lattice result for p(X) is 
displayed in Fig. |] for f3 = 2.0, V = 10 4 and f3 = 2.4, 
V = 16 4 , respectively. Note the strong decrease in (ipip) 
(in lattice units) for {3 = 2.4, cf. Eq. M) and Table §. 
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FIG. 1. Spectral density of the lattice Dirac operator for 
f3 — 2.0 and 2.4. Only positive eigenvalues are plotted. 

We are particularly interested in the region of small 
eigenvalues to check the predictions from chiral RMT. 
According to Ref. (Tsl, staggered fermions in SU(2) have 
the symmetries of the chGSE. Analytical expressions can 
be obtained in the framework of RMT for the micro- 
scopic spectral density and the distribution of the small- 
est eigenvalue by slight modifications of results computed 
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FIG. 2. Microscopic spectral density (upper row) and distribution of the smallest eigenvalue of the Dirac operator for different 
lattice parameters. From left to right the values of (5 are: 2.0, 2.0, 2.2, and 2.4. The histograms represent lattice data, the 
dashed curves are predictions from chiral RMT with Nf = v = 0. 



for Laguerre symplectic ensembles pp| , pq ] . Incorporating 
the chiral structure of the Dirac operator, we obtain from 
Ref. M 



»i r \ 

,2 / J„,„,2 



Ps(z) — 2z / duu / dv[J4 a -i(2uvz) J± a (2uz) 
Jo Jo 

-vJ ia -i(2uz)J ia (2uvz)] (3) 

with 4a = N t-\-2v+\, where Nt is the number of massless 
flavors and v is the topological charge. For our quenched 
data, 4a — 1 since v = as explained in the introduc- 
tion. According to Eq. (^), lattice data for p s {z) are 
constructed from the numerical eigenvalue density using 
a scale V(-ipip). This scale is determined by the data, 
hence the random-matrix predictions are parameter-free. 
Similarly, the distribution of the smallest eigenvalue for 
N f = v = follows from Ref. |S, 



-P(Amin) 



c(cA min ) 3/2 / 3/2 (cA min )e 



-i(cA mI „) 2 



(4) 



where c = Vfyij)) is the same scale as above. In Fig. g 
we have plotted the lattice results for p s (z) and P(A m i n ) 
together with the analytical results of Eqs. (|^) and (^) 
for four different combinations of (3 and lattice size. The 
agreement between lattice data and the parameter-free 
RMT predictions is impressive. Note that the RMT re- 
sults were derived in the limit V —* oo. Clearly, the 
agreement improves as the physical volume increases, i.e., 
with larger lattice size and smaller f3. From the results 
for (3 = 2.0 we observe that the agreement with RMT 
improves with increasing lattice size while the value of 



condensate remains the same. This suggests that a sim- 
ilar improvement will occur for (3 — 2.2 and (3 = 2.4. 
These values are just below the (3- value above which (ipip) 
approaches zero, where the above RMT results are inap- 
plicable, and an increased sensitivity to the size of the 
lattice is expected. Since P(A m i n ) for these couplings 
agrees with the RMT distribution for zero topological 
charge we expect that the discrepancy for p s {z) is not 
due to a superposition of configurations with different 
topological charge. We hope that future work will clarify 
this issue. 

Related quantities testing similar properties are the 
higher-order spectral correlation functions, in particular 
the two-point function which enters in the computation 
of scalar susceptibilities. The n-point correlation func- 
tion R n {x\, . . . , x n ) is defined as the probability density 
of finding a level (regardless of labeling) around each 
of the points x\, . . . , x n . The two-level cluster function 
T2(x, y), which contains only the non-trivial correlations, 
is defined by T2(x, y) = —R2(x, y) + Ri(x)R\ (y), i.e., the 
disconnected part is subtracted. For the chGUE, there 
are analytical arguments |27]] that the microscopic cor- 
relations are universal, and the same is expected for the 
chGSE. In this case, the predictions from RMT can again 
be obtained from the results of Ref. |^(| , but we do not 
write down the explicit expressions here. In Fig. ^|, we 
have plotted data for p s (x, y) for (3 = 2.0 on an 8 4 lattice 
as a function of x for some fixed value of y, along with the 
analytical random-matrix prediction. Clearly, the statis- 
tics are not as good as for the one-point function, but the 
agreement is still quite impressive. 

Finally, we have checked the Leutwyler-Smilga sum 
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to the continuum limit, and we expect random-matrix 
results to be a useful tool in the analysis. 
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a CRAY T90 at the Forschungszcntrum Julich and on a 
CRAY T3E at the HWW Stuttgart. 



FIG. 3. Microscopic limit of the two-level cluster function 
for some fixed value of y. The histogram represents data 
for (3 = 2.0 on an 8 4 lattice, the dashed curve shows the 
random-matrix prediction. 



rulc (E„ K 2 )/ y2 = W)7 2 appropriate for the chGSE 
Hps) . The numerical results are compared with the an- 
alytical predictions in Table |l| Again, the agreement im- 
proves with physical volume. 

TABLE I. Chiral condensate and a comparison of lattice 
data and analytical predictions for the Leutwyler-Smilga sum 
rule for \Z 2 . 







L 



(W>) 



<E„0/f 2 



W>72 



2.0 
2.0 
2.2 
2.4 



8 
10 
8 
16 



0.1228(25) 
0.1247(22) 
0.0556(19) 
0.00863(48) 



8.20(20) -lO -3 
7.97(30)-10~ 3 
1.67(03)-10~ 3 
3.97(14) -10~ 5 



7.54(31) -10 -3 
7.78(27) lO" 3 
1.55(11)-10 -3 
3.72(42) -10 -5 



In summary, we have performed a high-statistics study 
of the eigenvalue spectrum of the lattice QCD Dirac oper- 
ator with particular emphasis on the low-lying eigenval- 
ues. In the absence of a formal proof, our results provide 
very strong and direct evidence for the universality of p s . 
In the strong coupling domain, the agreement with an- 
alytical predictions from random matrix theory is very 
good. On the scale of the smallest eigenvalues, agree- 
ment is found even in the weak-coupling regime. Fur- 
thermore, we found that the microscopic two-level cluster 
function agrees nicely with random-matrix predictions 
and that the Leutwyler-Smilga sum rule for A~ 2 is sat- 
isfied more accurately with increasing physical volume. 
We predict that corresponding lattice data for SU(2) 
with Wilson fermions and for SU(3) with staggered and 
Wilson fermions will be described by random matrix re- 
sults for the GOE, chGUE, and GUE, respectively @. 
(The Ua(1) symmetry is absent for the Hermitean Wil- 
son Dirac operator.) The identification of universal fea- 
tures in lattice data is both of conceptual interest and of 
practical use. In particular, the availability of analytical 
results allows for reliable extrapolations to the chiral and 
thermodynamic limits. In future work we hope to ana- 
lyze the fate of the fermionic zero modes in the approach 
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